A PRIMAL BARVINOK ALGORITHM 
BASED ON IRRATIONAL DECOMPOSITIONS 



MATTHIAS KOPPE 

Abstract. We introduce variants of Barvinok's algorithm for counting 
lattice points in polyhedra. The new algorithms are based on irrational 
signed decomposition in the primal space and the construction of ra- 
tional generating functions for cones with low index. We give compu- 
tational results that show that the new algorithms are faster than the 
existing algorithms by a large factor. 



1. Introduction 

Twelve years have passed since Alexander Barvinok's amazing algo- 
rithm for counting lattice points in polyhedra was published (Barvinok, 
1994). In the mean time, efficient implementations (De Loera et al., 2004b, 
Verdoolaege et al., 2005) were designed, which helped to make Barvinok's al- 
gorithm a practical tool in many applications in discrete mathematics. The 
implications of Barvinok's technique, of course, reach far beyond the domain 
of combinatorial counting problems: For example, De Loera et al. (2005b) 
pointed out applications in Integer Linear Programming, and De Loera et al. 
(2006a, b) obtained a fully polynomial-time approximation scheme (FPTAS) 
for optimizing arbitrary polynomial functions over the mixed-integer points 
in polytopes of fixed dimension. 

Barvinok's algorithm first triangulates the supporting cones of all vertices 
of a polytope, to obtain simplicial cones. Then, the simplicial cones are re- 
cursively decomposed into unimodular cones. It is essential that one uses 
signed decompositions here; triangulating these cones is not good enough 
to give a polynomiality result. The rational generating functions of the 
resulting unimodular cones can then be written down easily. Adding and 
subtracting them according to the inclusion-exclusion principle and the the- 
orem of Brion (1988) gives the rational generating function of the polytope. 
The number of lattice points in the polytope can finally be obtained by 
applying residue techniques on the rational generating function. 

The algorithm in the original paper (Barvinok, 1994) worked explicitly 
with all the lower-dimensional cones that arise from the intersecting faces of 
the subcones in an inclusion-exclusion formula. Later it was pointed out that 
it is possible to simplify the algorithm by computing with full-dimensional 
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cones only, by making use of Brion's "polarization trick" (see Barvinok and 
Pommersheim, 1999, Remark 4.3): The computations with rational gener- 
ating functions are invariant with respect to the contribution of non-pointed 
cones (cones containing a non-trivial linear subspace). By operating in the 
dual space, i.e., by computing with the polars of all cones, lower-dimensional 
cones can be safely discarded, because this is equivalent to discarding non- 
pointed cones in the primal space. The practical implementations also rely 
heavily on this polarization trick. 

In practical implementations of Barvinok's algorithm, one observes that 
in the hierarchy of cone decompositions, the index of the decomposed cones 
quickly descends from large numbers to fairly low numbers. The "last mile," 
i.e., decomposing many cones with fairly low index, creates a huge number 
of unimodular cones and thus is the bottleneck of the whole computation in 
many instances. 

The idea of this paper is to stop the decomposition when the index of 
a cone is small enough, and to compute with generating functions for the 
integer points in cones of small index rather than unimodular cones. When 
we try to implement this simple idea in Barvinok's algorithm, as outlined 
in section 3, we face a major difficulty, however: Polarizing back a cone of 
small index can create a cone of very large index, because determinants of 
d X d matrices are homogeneous of order d. 

To address this difficulty, we avoid polarization altogether and perform 
the signed decomposition in the primal space instead. To avoid having to 
deal with all the lower-dimensional subcones, we use the concept of irrational 
decompositions of rational polyhedra. Beck and Sottile (2005) introduced 
this notion to give astonishingly simple proofs for three theorems of Stanley 
on generating functions for the integer points in rational polyhedral cones. 
Using the same technique. Beck et al. (2005) gave simplified proofs of the- 
orems of Brion and Lawrence-Varchenko. An irrational decomposition of 
a polyhedron is a decomposition into polyhedra whose proper faces do not 
contain any lattice points. Counting formulas for lattice points based on irra- 
tional decompositions therefore do not need to take any inclusion-exclusion 
principle into account. 

We give an explicit construction of a uniform irrational shifting vector s 
for a cone v + K with apex v such that the shifted cone {v + s) + K has 
the same lattice points and contains no lattice points on its proper faces 
(section 4). More strongly, we prove that all cones appearing in the signed 
decompositions of {v + s) + K in Barvinok's algorithm contain no lattice 
points on their proper faces. Therefore, discarding lower-dimensional cones 
is safe. Despite its name, the vector s only has rational coordinates, so after 
shifting the cone by s, large parts of existing implementations of Barvinok's 
algorithm can be reused to compute the irrational primal decompositions. 

In section 5, we show the precise algorithm. We also show that the same 
technique can be applied to the "homogenized version" of Barvinok's algo- 
rithm that was proposed by De Loera et al. (2004a). 

In section 6, we extend the irrationalization technique to non-simplicial 
cones. This gives rise to an "all-primal" Barvinok algorithm, where also 
triangulation of non-simplicial cones is performed in the primal space. This 
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allows us to handle problems where the triangulation of the dual cones is 
hard, e.g., in the case of cross polytopes. 

Finally, in section 7, we report on computational results. Results on 
benchmark problems show that the new algorithms are faster than the exist- 
ing algorithms by orders of magnitude. We also include results for problems 
that could not previously be solved with Barvinok techniques. 

2. Barvinok's algorithm 

Let P C R'' be a rational polyhedron. The generating function of P fl Z*^ 
is defined as the formal Laurent series 

5p(z) = z" e Z[[zi,...,Zrf,2;f\...,z^^]], 

ctGPnZ'* 

using the multi-exponent notation z" = nf=i-^r*- ^ bounded, gp is 
a Laurent polynomial, which we consider as a rational function gp. If P is 
not bounded but is pointed (i.e., P does not contain a straight line), there is 
a non-empty open subset U C such that the series converges absolutely 
and uniformly on every compact subset U to a rational function gp. If 
P contains a straight line, we set gp = 0. The rational function gp £ 
Q(zi, . . . , Zd) defined in this way is called the rational generating function 
of P n Z'^. 

Barvinok's algorithm computes the rational generating function of a poly- 
hedron P. It proceeds as follows. By the theorem of Brion (1988), the 
rational generating function of a polyhedron can be expressed as the sum 
of the rational generating functions of the supporting cones of its vertices. 
Let Vj G Q"' be a vertex of the polyhedron P. Then the supporting cone 
Vj + Ci of Vj is the (shifted) polyhedral cone defined by Vj + cone(P — Vj). 
Every supporting cone Vj + Ci can be triangulated to obtain simplicial cones 
Vj + Cij. Let K = V + BH'^ be a simplicial full-dimensional cone, whose 
basis vectors hi, . . . ,hd (i.e., representatives of its extreme rays) are given 
by the columns of some matrix B S z'^^'^. We assume that the basis vectors 
are primitive vectors of the standard lattice Z"^. Then the index of K is 
defined to be indK = \det B\; it can also be interpreted as the cardinality 
of nn Z"^, where 11 is the fundamental parallelepiped of K, i.e., the half-open 
parallelepiped 

n = v + {^ti Aib, :0< A, <l}. 

We remark that the set 11 n Z'^ can also be seen as a set of representatives 
of the cosets of the lattice BZ'^ in the standard lattice Z"^; we shall make 
use of this interpretation in section 3. Barvinok's algorithm now computes 
a signed decomposition of the simplicial cone K to produce other simplicial 
cones with smaller index. To this end, the algorithm constructs a vector 
w £ Z'^ such that 

w = aibi + • • • + ttdhd with \ai\ < |det Pj""^'''^ < 1. (1) 

This can be accomplished using integer programming or lattice basis reduc- 
tion. The cone is then decomposed into cones spanned by d vectors from the 
set {bi, . . . , b^;, w}; each of the resulting cones then has an index bounded 
above by {indK)^'^~^^/'^. In general, these cones form a signed decomposition 
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Figure 1. A triangulation of the cone of index 5 gener- 
ated by and into the two cones spanned by {b^,w} 
and {b^,w}, having an index of 2 and 3, respectively. 
We have the inclusion-exclusion formula 5'cone{bi,b2}(^) = 
5'cone{bi,w}(z) + S'concjba.w} (z) " S'conciw} (z) ; here the one- 
dimensional cone spanned by w needed to be subtracted. 



of K (see Figure 2); if w lies inside K, they form a triangulation of K (see 
Figure 1). The resulting cones and their intersecting proper faces (arising 
in an inclusion-exclusion formula) are recursively processed, until unimodu- 
lar cones, i.e., cones of index 1 are obtained. Finally, for a unimodular cone 
V + , the rational generating function can be easily written down as 

(2) 



where a is the unique integer point in the fundamental parallelepiped of the 
cone. We summarize Barvinok's algorithm below. 

Algorithm 1 (Barvinok's original (primal) algorithm). 

Input: A polyhedron P C R*^ given by rational inequalities. 
Output: The rational generating function for P fl Z*^ in the form 



5'p(z) 



where ej € {±1}, aj G Z*^, and bj 



(3) 



1. Compute all vertices Vj and corresponding supporting cones Ci of P. 

2. Triangulate Ci into simplicial cones Cij, keeping track of all the inter- 
secting proper faces. 

3. Apply signed decomposition to the cones Vj + Cij to obtain unimodular 
cones Vj + Ciji, keeping track of all the intersecting proper faces. 

4. Compute the unique integer point a^ in the fundamental parallelepiped 
of every resulting cone Vj + Ciji . 

5. Write down the formula (3). 
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Figure 2. A signed decomposition of the cone of in- 
dex 5 generated by and b^ into the two unimodu- 
lar cones spanned by {b^,w} and {b^,w}. We have the 
inclusion-exclusion formula fi-conclbi.bajlz) = 9conc{bi,w}(z) - 

9'cone{b2,w} (^) + 5cone{w} 

(z). 



The recursive decomposition of cones defines a decomposition tree. Due 
to the descent of the indices in the signed decomposition procedure, the 
following estimate holds for its depth: 

Lemma 2 (Barvinok, 1994). Let BTV^ be a simplicial full- dimensional cone, 
whose basis is given by the columns of the matrix B G Z'^^'^. Let D = |det-B|. 
Then the depth of the decomposition tree is at most 



Because at each decomposition step at most 0(2^^) cones are created and 
the depth of the tree is doubly logarithmic in the index of the input cone, 
Barvinok could obtain a polynomiality result in fixed dimension: 

Theorem 3 (Barvinok, 1994). Letd be fixed. There exists a polynomial-time 
algorithm for computing the rational generating function of a polyhedron 
P C R*^ given by rational inequalities. 

Later the algorithm was improved by making use of Brion's "polariza- 
tion trick" (see Barvinok and Pommersheim, 1999, Remark 4.3): The com- 
putations with rational generating functions are invariant with respect to 
the contribution of non-pointed cones (cones containing a non-trivial linear 
subspace) . The reason is that the rational generating function of every non- 
pointed cone is zero. By operating in the dual space, i.e., by computing with 
the polars of all cones, lower-dimensional cones can be safely discarded, be- 
cause this is equivalent to discarding non-pointed cones in the primal space. 



kiD) 



1 + 



log2 log2 D 



(4) 
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Algorithm 4 (Dual Barvinok algorithm). 

Input: A polyhedron P C R'^ given by rational inequalities. 
Output: The rational generating function for P PI Z*^ in the form 



(5) 



where G {±1}, e Z"^, and bjj G Z'^. 

1. Compute all vertices Vj and corresponding supporting cones Cj of P. 

2. Polarize the supporting cones Cj to obtain C°. 

3. Triangulate C° into simplicial cones C^?-, discarding lower-dimensional 
cones. 

4. Apply Barvinok's signed decomposition to the cones Vj + C°j to obtain 
cones Vj + C°ji, stopping decomposition when a unimodular cone is ob- 
tained. Discard all lower-dimensional cones. 

5. Polarize back C°ji to obtain cones Ciji. 

6. Compute the unique integer point in the fundamental parallelepiped 
of every resulting cone Vj + Ciji. 

7. Write down the formula (5). 

This variant of the algorithm is much faster than the original algorithm 
because in each step of the signed decomposition at most d, rather than 
0(2^^), cones are created. The practical implementations LattE (De Loera 
et al., 2004b) and barvinok (Verdoolaege et al., 2005) also rely heavily on 
this polarization trick. 

3. The Barvinok algorithm with stopped decomposition 

We start out by introducing a first variant of Barvinok's algorithm that 
stops decomposing cones before unimodular cones are reached. As we will 
see in the computational results in section 7, already the simple modification 
that we propose can give a significant improvement of the running time for 
some problems, at least in low dimension. 

Algorithm 5 (Dual Barvinok algorithm with stopped decomposition). 

Input: A polyhedron P C R'^ given by rational inequalities; the maximum 
index i. 

Output: The rational generating function for P H Z'^ in the form 



where € {±1}, Ai C Z'^ with \Ai\ < £, and hij G Z'^. 

1. Compute all vertices Vj and corresponding supporting cones Cj of P. 

2. Polarize the supporting cones Cj to obtain C°. 

3. Triangulate C° into simplicial cones C°j, discarding lower-dimensional 
cones. 

4. Apply Barvinok's signed decomposition to the cones Vj + C°j to ob- 
tain cones Vj + C°ji, stopping decomposition when a polarized-back cone 
Ciji = has index at most i. Discard all lower-dimensional cones. 

5. Polarize back C?-, to obtain cones Cj,;. 
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6. Enumerate the integer points in the fundamental parallelepipeds of all 
resulting cones Vj + Ciji to obtain the sets Ai. 

7. Write down the formula (6). 

As mentioned above, the integer points in the fundamental parallelepiped 
of a cone Vj + BijiR!^ can be interpreted as representatives of the cosets 
of the lattice BijiT/^ in the standard lattice Z'^. Hence they can be easily 
enumerated in step 6 by computing the Smith normal form of the generator 
matrix Biji] see Lemma 5.2 of Barvinok (1993). The Smith normal form can 
be computed in polynomial time, even if the dimension is not fixed (Kannan 
and Bachem, 1979). 

We remark that both triangulation and signed decomposition are done 
in the dual space, but the stopping criterion is the index of the polarized- 
back cones (in the primal space). The reason for this stopping criterion is 
that we wish to control the maximum number of points in the fundamental 
parallelepipeds that need to be enumerated. Indeed, when the maximum £ 
is chosen as a constant or polynomially in the input size, then Algorithm 5 
clearly runs in polynomial time (in fixed dimension). 

Each step of Barvinok's signed decomposition reduces the index of the 
decomposed cones. When the index of a cone 0°^^ is A, in the worst case 

the polarized-back cone Ciji has index A*^"^, where d is the dimension. If 
the dimension is too large, the algorithm often needs to decompose cones 
down to a very low index or even index 1, so the speed-up of the algorithm 
will be very limited. This can be seen from the computational results in 
section 7. 

4. Construction of a uniform irrational shifting vector 

In this section, we will give an explicit construction of an irrational shifting 
vector s for a simplicial cone w + K with apex v such that the shifted cone 
(v + s) + IT has the same lattice points and contains no lattice points on 
its proper faces. The "irrationalization" (or perturbation) will be uniform 
in the sense that also every cone arising during the Barvinok decomposition 
does not contain any lattice points on its proper faces. This will enable us 
to perform the Barvinok decomposition in the primal space, discarding all 
lower-dimensional cones. 

To accomplish this goal, we shall first describe a subset of the stability 
region of a cone ^^ + K with apex at v, i.e., the set of apex points v such 
that ^r + K contains the same lattice points as v + see Figure 3. 

Lemma 6 (Stability cube). Let v + WR!^ be a simplicial full- dimensional 
cone with apex at^i ^ Q*^, whose basis is given by the columns of the matrix 
B £ z*^^*^. Let h^, . . . ,b^ be a basis of the dual cone, given by the columns 
of the matrix B* = —{B~^)'^. 

Let D = |deti?|. Let A G Q"' and X £ Q'^ be defined by 

Ai = (b*,v) and = ^[DA^J + ^ for i = 1, . . . ,d. 

Let 

1 

V = —BX and o = — ■ — . 
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Figure 3. The stability region of a cone 



Then, for every v with ||v — v||^ < p, the cone v + BH'^ contains the same 
integer points as the cone v + BH'^ and does not have integer points on its 
proper faces. 

In the proof of the lemma, we will use of the inequality description (H- 
representation) of the simplicial cone v + BH'^. It is given by the basis 
vectors of the dual cone: 

V + BK'i = I X € R*^ : (b*, x) < (b*, v) for i = 1, . . . , } . (7) 
Proof of Lemma 6. Let A be defined by Aj = (b|, v). Then we have 

|A. - A,| < IIKIIi • ||v - v|U < IIKIIi • P < ^- (8) 

By (7), a point x € Z'^ lies in the cone v + B'R!^ if and only if 
(b*, x) < (b*, v) = Ai for i = 1, . . . , d. 

Likewise, x € v + BYi!^ if and only if 

(b*,x) < (b*,v) = Ai fori = l,...,d. 

Note that for x € Z*^, the left-hand sides of both inequalities are an integer 
multiple of -jj. Therefore, we obtain equivalent statements by rounding down 
the right-hand sides to integer multiples of . For the right-hand side of (4) 
we have by (8) 

A, = A, + (A, - A.) < (l^A,J + 1^ + ^ = 1 [D\\ + 1, (9a) 
\i = \ + (A, - A.) > (l^A,J + 1^ - ^ = 1 [DK\ , (9b) 

so Aj and Aj are rounded down to the same value [i^A^ J . Thus, the cone 
V -|- BYii^ contains the same integer points as the cone v + i?R^|_. Moreover, 
since the inequalities (9) are strict, the cone v -|- B'Ri^ does not have integer 
points on its proper faces. □ 
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For non-simplicial cones, we will give an algorithmic construction for a 
stability cube in section 6. 

Next we make use of the estimate for the depth of the decomposition tree 
in Barvinok's algorithm given in Lemma 2. On each level of the decom- 
position, the entries in the basis matrices can grow, but not by much. We 
obtain: 

Lemma 7. Let BH'^ be a simplicial full- dimensional cone, whose basis is 
given by the columns of the matrix B G 7/^^'^. Let D = |det-B|. Let C E Z_|_ 
he a number such that \Bi^j\ < C . 

Then all the basis matrices B of the cones that appear in the recursive 
signed decomposition procedure of Barvinok 's algorithm applied to BIV^ have 
entries bounded above by d^^^^C , where k{D) is defined by (4). 

Proof. Given a cone spanned by the columns bi, . . . ,b(i of the matrix B, 
Barvinok's algorithm constructs a vector w € Z'^ such that 

w = aibi H h aahd with \ai\ < |det < 1. (10) 

Thus |Iw||^ < dC. The cone is then decomposed into cones spanned by d 
vectors from the set {bi, . . . , b^, w}. Thus the entries in the corresponding 
basis matrices are bounded by dC. The result follows then by Lemma 2. □ 

If we can bound the entries of an integer matrix with non-zero determi- 
nant, we can also bound the entries of its inverse. 

Lemma 8. Let B € Z^^^ be a matrix with \Bij\ < C. Let D = jdeti?|. 
Then the absolute values of the entries of B~^ are bounded above by 

^{d-l)\C'~\ 

Proof. We have |(i?^"'^)fc,«| = |deti?(fc ;)|, where -B(fc,i) is the matrix ob- 
tained from deleting the k-th. row and Z-th column from B. Now the desired 
estimate follows from a formula for det B(^i^ i-j and from \Bij\ < C. □ 

Thus, we obtain a bound on the norm of the basis vectors of the polars 
of all cones occurring in the signed decomposition procedure of Barvinok's 
algorithm. 

Corollary 9 (A bound on the dual basis vectors). Let BIV^ be a simplicial 
full- dimensional cone, whose basis is given by the columns of the matrix 
B G Z'^^'^. Let D = |detS|. Let C be a number such that \Bij\ < C. 

Let B* = —{B~^)^ be the basis matrix of the polar of an arbitrary cone 
BYl% that appears in the recursive signed decomposition procedure applied 
to i?R^. Then, for every column vector b* of B* we have the estimate 

\\deiB- h*\\^ <{d- 1)1 (^d'^^'^Cy'' =: L (11) 
where k{D) is defined by (4). 

Proof. By Lemma 7, the entries of B are bounded above by d^^^'^C. Then 

the result follows from Lemma 8. □ 

The construction of the "irrational" shifting vector is based on the fol- 
lowing lemma. 



10 



MATTHIAS KOPPE 



Lemma 10 (The irrational lemma). Let M G Z_|_ he an integer. Let 

(2M'pff (2Mp) ■ ^^^^ 

Then (c,q) ^ Z for every c G Z'^ \ {0} with \\c\\^ < M. 

Proof. Follows from the principle of representations of rational numbers in 
a positional system of base 2M. □ 

Theorem 11. Let^ + BYHi^ he a simplicial full- dimensional cone with apex 
at w ^ Q'^, whose basis is given by the columns of the matrix B G Z"'^'^. 
Let D = |deti?|, C he a number such that < C and let v G Q'^ and 

p G Q+ be the data from Lemma 6 describing the stability cube o/v + 5R!|. 
Let < r G Z such that r^^ < p. Using 



^ log2 log2 D 
log2 " 



(13) 



(2M)i' (2M)2''"' (2M)" 



'2 

L = {d- iy.{d''Cf-^, and M = 2L, define 

1/1 1 
s = - 

r 

Finally let v = v + s. 

(i) VKe /laue (v + ^R^:^) n Z"' = (v + BR'l) n Z"^, i.e., the shifted cone has 
the same set of integer points as the original cone. 

(ii) The shifted cone v + i3R5| contains no lattice points on its proper faces. 

(iii) More strongly, all cones appearing in the signed decompositions of the 
shifted cone v + BH'^ in Barvinok 's algorithm contain no lattice points 
on their proper faces. 

Proof. Part (i). This follows from Lemma 6 because v clearly lies in the 
open stability cube. 

Parts (ii) and (iii). Every cone appearing in the course of Barvinok's 
signed decomposition algorithm has the same apex v as the input cone and 
a basis B G Z'^^'^ with |det-B| < D. Let such a S be fixed and denote by 
h* the columns of the dual basis matrix B* = —{B~^)~^. Let z G Z^ be an 
arbitrary integer point. We shall show that z is not on any of the facets of 
the cone, i.e., 

(b*,z-v)7^0 foT i = l,...,d. (14) 
Let i G {1, . . . , c?} arbitrary. We will show (14) by proving that 

(detS -b^v) ^ Z. (15) 

Clearly, if (15) holds, we have (b*,v) ^ (deti?)~^Z. But since (b*,z) G 
(det -B)~"'^Z, we have (b*,z — v) ^ Z; in particular it is nonzero, which 
proves (14). 

To prove (15), let c = det^-b*. By Corollary 9, we have ||c||^ < L < M. 
Now Lemma 10 gives (c, s) ^ ^Z. Observing that by the definition (6), we 
have 

(c,v) = {c,-BX) G iZ. 

Therefore, we have (c,v) = (c,v + s) ^ ^Z. This proves (15), and thus 
completes the proof. □ 
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5. The irrational algorithms 
The following is our variant of the Barvinok algorithm. 
Algorithm 12 (Primal irrational Barvinok algorithm). 

Input: A polyhedron P C R"^ given by rational inequalities; the maximum 
index (.. 

Output: The rational generating function for P fl Z'^ in the form 

nj=i(i-z'''0 

where G {±1}, Ai C Z'^ with \Ai\ < i, and hij G Z'^. 

1. Compute all vertices Vj and corresponding supporting cones of P. 

2. Polarize the supporting cones Cj to obtain C°. 

3. Triangulate C° into simplicial cones C°j, discarding lower-dimensional 
cones. 

4. Polarize back C^?- to obtain simplicial cones Cij. 

5. Irrationalize all cones by computing new apex vectors Vjj G Q"^ from Vj 
and Cij as in Theorem 11. 

6. Apply Barvinok's signed decomposition to the cones Vij + Cij, discarding 
lower-dimensional cones, until all cones have index at most i. 

7. Enumerate the integer points in the fundamental parallelepipeds of all 
resulting cones to obtain the sets Ai. 

8. Write down the formula (16). 

Theorem 13. Algorithm 12 is correct and runs in polynomial time when the 
dimension d is fixed and the maximum index i is bounded by a polynomial 
in the input size. 

Proof. This is an immediate consequence of the analysis of Barvinok's algo- 
rithm. The irrationalization (step 5 of the algorithm) increases the encoding 
length of the apex vector only by a polynomial amount, because the dimen- 
sion d is fixed and the depth k only depends doubly logarithmic on the initial 
index of the cone. □ 

The same technique can also be applied to the "homogenized version" of 
Barvinok's algorithm that was proposed by De Loera et al. (2004a); see also 
De Loera et al. (2004b, Algorithm 11). 

Algorithm 14 (Irrational homogenized Barvinok algorithm). 

Input: A polyhedron P C R'^ given by rational inequalities in the form 
Ax < b; the maximum index i. 

Output: A rational generating function in the form (16) for the integer points 
in the homogenization of P, i.e., the cone 

C={(^x,0:xGP, eGR+}. (17) 

1. Consider the inequality description for C; it is given by Ax — b,^ < 0. 
The polar C° then has the rays {Ai^., —bi), i = 1, . . . ,m. 

2. Triangulate C° into simplicial cones Cj, discarding lower-dimensional 
cones. 

3. Polarize back the cones Cj to obtain simplicial cones Cj. 
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4. Irrationalize the cones Cj to obtain shifted cones Vj + Cj. 

5. Apply Barvinok's signed decomposition to the cones + Cj, discarding 
lower-dimensional cones, until all cones have index at most i. 

6. Write down the generating function. 

6. Extension to the non-simplicial case 

For polyhedral cones with few rays and many facets, it is usually much 
faster to perform triangulation in the primal space than in the dual space, 
cf. Biieler et al. (2000). In this section, we show how to perform both 
Barvinok decomposition and triangulation in the primal space. 

The key idea is to use linear programming to compute a subset of the 
stability region of the non-simplicial cones. 

Lemma 15. There is a polynomial-time algorithm that, given the vertex 
V e Q*^ and the facet vectors b* € Z'^, i = 1, . . . ,m, of a full- dimensional 
polyhedral cone C = v -|- i?R" , where n > d, computes a point v € Q*^ 
and a positive scalar p € Q such that for every v in the open cube with 
||v — v||(^ < p, the cone v -|- i?R" has no integer points on its proper faces 
and contains the same integer points as v -\- i?R" . 

Proof. We maximize p subject to the linear inequalities 

{hlv) + \\h*\\iP< L(b*,v)J +1, (18a) 
-(b*,v) + ||b:||ip<-L(K,v)J, (18b) 

where v € R'^ and p G R-)_. We can solve this linear optimization problem 
in polynomial time. Let (v, p) be an optimal solution. Let v G R"' with 
||v — v||^ < p. Let X G (v -|- BH'^) H Z"'. Then we have for every i G 
{1, . . . ,m} 

(b*,x)<(b*,v) 

= (b*,v) + (b*,v-v) 
<(b*,v) + ||b*|U|v-v||^ 

< (K,v) + ||b*||ip 

< L(b*,v)J+l by (18a). 

Because (b*,x) is integer, we actually have (b*,x) < [(b*, v)J. Thus, x lies 
in the cone v -|- BJH^. Conversely, let x G (v -|- BIV^) n Z'^. Then, for every 
i G {1, . . . , m}, we have (b*, x) < (b*, v). Since x G Z"', we can round down 
the right-hand side and obtain 

(K,x)< L(b*,v)j 

<(b*,v)-||b*||,p by (18b) 

< (b*,v) - ||b*||J|v-v||^ 
<(b,r,v) + (b:,v-v) 

= (b,r,v). 

Thus, X G V -|- BYV^. Moreover, since the inequality is strict, x does not lie 
on the face (b*, x) = (b*, v) of the cone v + i?R^[. □ 
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Lemma 16 (Bound on the index of all subcones). Let bj G Z'', i = 1, . . . ,n, 
be the generators of a full- dimensional polyhedral cone K CR'^. Then the 
cones of any triangulation of K have an index hounded by 

= fmax^^i llbifj . (19) 

Proof. Let B G Z'^^'^ be the generator matrix of a full-dimensional cone of 
a triangulation of K; then the columns B form a subset {bj^, . . . ,bi^} C 
{bi, . . . , b„}. Therefore 

Idet^l < Yl \\hij < (maxti ||bif)"^^ 
fc=i 

giving the desired bound. □ 

With these preparations, the following corollary is immediate. 

Corollary 17. Letv + BIV^ be a full- dimensional polyhedral cone with apex 
atvG Q"', whose basis is given by the columns of the matrix B G Z'^^"'. Let 
D be defined by (19). Let v G Q"^ and p G Q+ be the data from Lemma 15 
describing the stability cube o/v + i?R". Using these data, construct v as 
in Theorem IL Then the assertions of Theorem, 11 hold. 

Algorithm 18 (All-primal irrational Barvinok algorithm). 

Input: A polyhedron P C R*^ given by rational inequalities; the maximum 
index (.. 

Output: The rational generating function for P n Z'^ in the form (16). 

1. Compute all vertices Vj and corresponding supporting cones Ci of P. 

2. Irrationalize all cones by computing new apex vectors Vj G Q'^ from Vj 
by Corollary 17. 

3. Triangulate Vj + Cj into simplicial cones Vj + Cjj, discarding lower- 
dimensional cones. 

4. Apply Barvinok's signed decomposition to the cones Vj + Cij, until all 
cones have index at most L 

5. Enumerate the integer points in the fundamental parallelepipeds of all 
resulting cones to obtain the sets Ai. 

6. Write down the formula (16). 

7. Computational experiments 

Algorithms 12 and 18 have been implemented in a new version of the 
software package LattE, derived from the official LattE release 1.2 (De Loera 
et al., 2005a). The new version, called LattE macchiato, is freely available on 
the Internet (Koppe, 2006). In this section, we discuss some implementation 
details and show the results of first computational experiments. 

7.1. Two substitution methods. When the generating function gp has 
been computed, the number of lattice points can be obtained by evaluating 
gp(l). However, 1 is a pole of every summand of the expression 
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The method implemented in LattE 1.2 (see De Loera et al., 2004b) is to 
use the polynomial substitution 



for a suitable vector A. Then the constant coefficient of the Laurent expan- 
sion of every summand about s = is computed using polynomial division. 
The sum of all the constant coefficients finally gives the number of lattice 
points. 

Another method from the literature (see, for instance, Barvinok and Pom- 
mersheim, 1999) is to use the exponential substitution 



for a suitable vector A. By letting r — > 0, one then obtains the formula 



where td^-fc is the so-called Todd polynomial. In LattE macchiato, the 
exponential substitution method has been implemented in addition to the 
existing polynomial substitution; see De Loera and Koppe (2006) for imple- 
mentation details. 

7.2. Implementation details. We enumerate the lattice points in the fun- 
damental parallelepiped by computing the Smith normal form of the gener- 
ator matrix B; see Lemma 5.2 of Barvinok (1993).^ For computing Smith 
normal forms, we use the library LiDIA, version 2.2.0. For solving the linear 
program in Lemma 15, we use the implementation of the revised dual sim- 
plex method in exact rational arithmetic in cddlib, version 0.94a (Fukuda, 
2005). All other computations are done using the libraries NTL, version 5.4 
(Shoup, 2005) and GMP, version 4.1.4 for providing exact integer and ra- 
tional arithmetic. 

7.3. Evaluation of variants of the algorithms. We compare the vari- 
ants of the algorithms using test instances that can also be solved without 
the proposed irrationalization techniques. We consider the test instances 
hickerson-12, hickerson-13, and hickerson-14, related to the manu- 
script by Hickerson (1991). They describe simplices in and that 
contain 38, 14, and 32 integer points, respectively. The examples are good 
test cases for our algorithms because the vertices and cones are trivially 
computed, and all computation time is spent in the Barvinok decomposi- 
tion. We show the results in Table 1, Table 2, and Table 3. The tables show 
results for the following methods: 

1. Methods without irrationalization, using polarization to avoid comput- 
ing with lower-dimensional cones: 

(a) LattE 1.2 (De Loera et al., 2005a), decomposing down to unimodular 
cones in the dual space (Algorithm 5 with £ = 1). 

(b) Likewise, but using the implementation in the library barvinok by 
Verdoolaege (2006), version 0.21. 

-'^The author wishes to thank Susan Marguhes for prototyping the enumeration code. 



z = 



z = 



(exp{rAi}, . . . ,exp{TArf}) 




(20) 
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Without irrationalization With irrationahzation 



Time (s) Time (s) 



Max. 




LattE 


barv . 


LattE macchiato 




LattE macchiato 


index 


Cones 


V 1.2 


v0.21 


Poly 


Exp 


Cones 


Poly 


Exp 


1 


11625 


17.9 


11.9 


10.0 


16.7 


7929 


7.8 


12.7 


10 


4251 






6.9 


7.0 


803 


1.9 


1.6 


100 


980 






6.9 


2.1 


84 


1.3 


0.3 


200 


550 






9.1 


1.5 


76 


1.3 


0.3 


300 


474 






9.9 


1.4 


58 


1.4 


0.3 


500 


410 






11.7 


1.3 


42 


1.6 


0.3 


1000 


130 






7.2 


0.7 


22 


1.7 


0.2 


2000 


7 






2.2 


0.2 


22 


1.8 


0.2 


5000 


7 






2.8 


0.2 


7 


2.8 


0.2 



(c) LattE macchiato, decomposing cones in the dual space, until all 
cones in the primal space have at most index i (Algorithm 5), then 
using polynomial substitution. We show the results for different 
values of I. 

(d) Likewise, but using exponential substitution. 

2. Methods with irrationalization, performing triangulation in the dual 
space and Barvinok decomposition in the primal space (Algorithm 12): 

(a) LattE macchiato with polynomial substitution. 

(b) LattE macchiato with exponential substitution. 

The table shows computation times in CPU seconds on a PC with a Pen- 
tium M processor with 1.4 GHz. It also shows the total number of simplicial 
cones created in the decomposition, using the different variants of LattE; 
note that we did not measure the number of simplicial cones that the li- 
brary barvinok produced. 

We can make the following observations: 

i. By stopping Barvinok decomposition before the cones are unimodular, 
it is possible to significantly reduce the number of simplicial cones. 
This effect is much stronger with irrational decomposition in the primal 
space than with decomposition in the dual space. 

ii. The newly implemented exponential substitution has a computational 
overhead compared to the polynomial substitution that was imple- 
mented in LattE 1.2. 

iii. However, when we compute with simplicial, non-unimodular cones, the 
exponential substitution becomes much more efficient than the poly- 
nomial substitution. Hence the break-even point between enumeration 
and decomposition is reached at a larger cone index. 

The reason is that the inner loops are shorter for the exponential sub- 
stitution; essentially, only a sum of powers of scalar products needs to 
be evaluated in the formula (20). This can be done very efficiently. 

iv. The best results are obtained with the irrational primal decomposition 
down to an index of about 500 to 1000 and exponential substitution. 
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Table 2. Results for hickerson-13 



Without irrationalization With irrationahzation 









Time (s) 






Time (s) 


iVldA.. 






barv . 


LattE macchiato 




LattE macchiato 


indsx 


1 v~vn eta 


V 1 . Z 


v0.21 


Poly 


Exp 


v^oncs 


Poly 


Exp 


1 


465 540 


795 


589 


421 


707 


483 507 


479 


770 


10 


272 922 






345 


428 


55 643 


117 


109 


100 


142 905 






489 


249 


9158 


83 


22 


200 


122 647 






625 


222 


6150 


93 


17 


300 


98 654 






903 


199 


4674 


105 


14 


500 


90 888 






1056 


193 


3 381 


137 


13 


1000 


73 970 






1648 


190 


2 490 


174 


13 


2000 


66 954 






2166 


201 


1857 


237 


14 


5000 


49168 






5040 


286 


1488 


354 


18 


10000 


43511 






7278 


370 


1011 


772 


34 



Table 3. Results for hickerson-14 



Without irrationalization With irrationalization 

Time (s) Time (s) 



LattE macchiato LattE macchiato 
barv . 



index 


Cones vl.2 


v0.21 Poly 


Exp 


Cones 


Poly 


Exp 


1 


1682 743 4017 


15 284 2 053 


3466 


552 065 


792 


1244 


10 


1027619 


1736 


2177 


49 632 


168 


143 


100 


455 474 


2 294 


1089 


8 470 


128 


29 


200 


406491 


2 791 


990 


5 554 


157 


22 


300 


328340 


4131 


875 


4 332 


187 


19 


500 


303 566 


4911 


842 


3464 


235 


18 


1000 


236 626 


8 229 


807 


2 384 


337 


18 


2000 


195 368 


12122 


817 


1792 


481 


21 


5000 


157496 


22 972 


1034 


1276 


723 


27 


10000 


128 372 


31585 


1270 


956 


1095 


38 



7.4. Results for challenge problems. In Table 4 we show the results for 
some larger test cases related to Hickerson (1991). We compare LattE 1.2 
with our implementation of irrational primal decomposition (Algorithm 12) 
with maximum index 500. The computation times are given in CPU sec- 
onds. The computations with LattE 1.2 were done on a PC Pentium M, 
1.4 GHz; the computations with LattE macchiato were done on a slightly 
slower machine, a Sun Fire V890 with UltraSPARC-IV processors, 1.2 GHz. 

Both the traditional Barvinok algorithm (Algorithm 5 with i = 1) and 
the homogenized variant of Barvinok's algorithm (De Loera et al., 2004a) 
do not work well for cross polytopes. The reason is that triangulation is 
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Table 4. Results for larger Hickerson problems 



n 


d 


Lattice 
points 


LattE vl.2 


LattE macchiato 


Cones 


Time 


Cones 


Time 


15 


7 


20 


293 000 


10 min 55 s 


2 000 


22 s 


16 


8 


54 


3 922 000 


3 h 35 min 


19 000 


3 min 56 s 


17 


8 


18 






2 655 000 


7 h 59 min 


18 


9 


44 


61500 000 


77 h 00 min 


200 000 


49min 12s 


20 


10 


74 






2 742 000 


13 h 05 min 



Table 5. Results for cross polytopes 



d 


Without irrationalization 


All-primal 


irrational 


Cones 


Time (s) 


Cones 


Time (s) 


4 


384 


1.1 




0.9 


5 


3 840 


6.5 




1.4 


6 


46 264 


91.7 




2.7 


7 


653 824 


1688.7 




5.5 


8 






1000 


12.3 


9 






2 000 


29.6 


10 






5 000 


74.8 


11 






11000 


189.1 


12 






24000 


483.0 


13 






53 000 


1231.2 


14 






114 000 


3145.6 


15 






245 000 


8180.9 



done in the dual space, so hypercubes need to be triangulated. We show the 
performance of the traditional Barvinok algorithm in Table 5. We also show 
computational results for the all-primal irrational algorithm (Algorithm 18 
with i = 500), using exponential substitution. The computation times are 
given in CPU seconds on a Sun Fire V440 with UltraSPARC-IIIi processors, 
1.6 GHz. 

A challenge problem related to the paper by Beck and Ho§ten (2006), 
case m = 42, could be solved using the all-primal irrational decomposition 
algorithm (Algorithm 18) with exponential substitution. The method de- 
composed the polyhedron to a total of 1.1 million simplicial cones of index 
at most 500. The computation took 66 000 CPU seconds on a Sun Fire 
V440 with UltraSPARC-IIIi processors, 1.6 GHz. The problem could not be 
solved previously because the traditional algorithms first tried to triangulate 
the polar cones, which does not finish within 17 days of computation. 

Conclusions and future work 

The above computational results with our preliminary implementation 
have shown that the proposed irrationalization techniques can speed up the 
Barvinok algorithm by large factors. 

A further speed-up can be expected from a refined implementation. For 
example, the choice of the irrational shifting vector is based on worst-case 
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estimates. It may be worthwhile to implement a randomized choice of the 
shifting vector (within the stability cube), using shorter rational numbers 
than those constructed in the paper. The randomized choice, of course, 
would not give the same guarantees as our deterministic construction. How- 
ever, it is easy and efficient to check, during the decomposition, if the gen- 
erated cones are all irrational; when they are not, one could choose a new 
random shifting vector (or resort to the one constructed in this paper) and 
restart the computation. 
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